"DSP chips can produce random numbers using proven algorithm"

Paul Mennen, Tektronix Inc
EDN January 21, 1991, page 141

Psuedorandom number sequences prove useful for many measurement techniques. The basic idea is to excite the system under test with a signal that simultaneously has all frequencies of interest instead of at one frequency at a time as in the swept-sine method. You then compare the Fourier transforms of the system input and output to generate a frequency response.

Of the hundreds of methods in use for generating psuedorandom sequences, this listing for the WE DSP16 use the Linear Congruential Sequence" (LCS). One advantage of the LCS is that it is simple enough that mathematicians have been able to prove many of its properties.

The defining equation for the LCS is:

R' = (aR + c) mod m

where R is the random number, a and c are constants, and mod m is used to signify the operation of taking the remainder after dividing by m (modulus). First pick a starting value for R, called the seed. (i.e. assign the seed as the first value for R). Then multiply R by a constant, add a different constant to the product, divide the sum by m, and assign the remainder resulting from this division to be the next value of R. These steps are repeated every time a new random number is required.

A little thought will reveal that the numbers produced (R) are between 0 and m-1 and that the sequence can not be longer than m. (assuming that a, c, m, and the seed are positive). Now you are only left with the non-trivial task of picking appropriate values for a, c, and m. (The choice of the seed turns out to be unimportant). I chose these numbers to satisfy the following 6 criteria. (Additionally define n = sqrt(m)).

1.) c and m are relatively prime

2.) a-1 is a multiple of every prime factor of m

3.) a-1 is a multiple of 4 if m is a multiple of 4

4.) a and c are both greater than n

5.) trunc(a/n) is a power of 2

6.) a mod n is less than n/2

The first three criteria insure that the sequence will be maximal length, ie the sequence will contain all the numbers between 0 and m-1 before repeating. (This result can be proved, but it's beyond the scope of the article.) Criterion 4 insures sufficient spectral flatness and distribution uniformity. Criteria 5 and 6 will turn out to simplify the computation without compromising the sequence quality.

Note that you can simplify the computation by choosing m to be a power of 2, since the division can be implemented by shifting. The DSP16 includes a 16-bit multiplier, so it would seem that m = 2**16 would be a logical choice. Assuming a maximal length sequence and a 512-kHz output rate, the sequence will repeat every .128 seconds (65536 / 512000). Too short for most practical applications.

The next logical choice then is m = 2**32, resulting in a sequence length of 2.33 hours (65536 * .128 seconds). This is long enough to meet even the most demanding situations.

For this choice of m, you need to do a 32-bit multiply. This operation requires multiple-precision arithmetic which complicates the code considerably. Define alo and ahi to be the lower and higher 16 bits of the constant a. Similarly denote the two halves of the random number as Rlo and Rhi. Now the LCS using multiple-precision arithmetic is

EQUATION 1:

R = c + alo*Rlo + ((alo*Rhi + ahi*Rlo + zz) << 16)

where

* signifies a product of two 16-bit numbers (32-bit result)
+ signifies a sum of two 32-bit numbers
<< signifies a logical left shift operation
zz = alo if Rlo is greater than 32768
zz = 0 otherwise

To see where this comes from first look at the multiple-precision product of the 32-bit constant "a" and the 32-bit random number "R". (Fig 5 in article).

Because the remainder after dividing by 2**32 (m) is simply the low order half of the 64 bit product, we don't even need to compute the upper half. Thus we need only save the lower half of the 2nd and 3rd partials and we don't need to compute the 4th partial at all. Because of criterion #5, the 3rd partial is a simple left shift, leaving only 2 required multiplies.

By now you should be able to verify equation 1 except for the mysterious zz term. This term stems from the fact that the partial products must be computed using unsigned multiplies for the 64-bit multiple-precision product to be correct. However, most multipliers (including the DSP16's) are 2's complement signed multipliers. Using a signed multiply for the 2nd, 3rd, and 4th partial products produces errors only in the unused high-order half of the product, so only the 1st partial remains a concern. Although alo will always be interpreted by the multiplier as a positive number (because of criterion #6), Rlo can be (and will be) any combination of 16 bits. Whenever the most significant bit of Rlo is 1, we need to correct the result of the signed multiply (with the zz term) to fool the multiplier into producing an unsigned result.

The actual values I chose for constants a and c are:

a = 107465 (hex) c = 234567 (hex)

You can verify that these constants satisfy the above 6 criteria. Note however that there is a large number of equally fine choices. Using these values for a and c equation 1 in psuedo code is

TEMP = 7465 * Rhi ; 2nd partial product (TEMP is 32 bits)
TEMP = TEMP + Rlo << 4
; add 3rd partial product (ahi = 10 hex) if (Rlo > 7fff)
; if Rlo > 32767 correct the 1st partial TEMP = TEMP + 7465
; product for using a signed multiply
TEMP = TEMP << 10 ; shift left 16 (discard 16 high order bits)
TEMP = TEMP + 7465 * Rlo ; add 1st partial product
R = TEMP + 234567 ; 32 bit add produces the new random number

All constants in the psuedo code are in hexidecimal. As before, Rlo and Rhi are the lower and upper 16 bits respectively of the random number R

For those familiar with the AT&T DSP16, the code segment which computes the LCS sequence follows. The computations have been carefully ordered to maximize the number of parallel operations. Note that the routine computes a 32-bit random number, but my application needed only 16 bits to send to a DAC. I chose to send the upper 16 bits the DAC since it turns out to be "more random" than the lower half.

Note that it takes only 15 cycles to output each random number. This means that if you were to use the fastest version of the DSP16A (25 ns) and the DSP chip was doing nothing else, you could output a new random number every 375 ns.

__________________________________________________________________________

           Code sequence for computing the LCS on the DSP16
__________________________________________________________________________
.ram
ALO:	int			/* storage for a(lo) - Multiplier	*/
CHI:	int			/* storage for c(hi) - Add constant	*/
CLO:	int			/* storage for c(lo)			*/
.endram
alo:	int 0x7465		/* a(lo)				*/
chi:	int 0x0023		/* c(hi)				*/
clo:	int 0x4567		/* c(lo)				*/
	i = 0			/* used as table increment		*/
	a0 = i			/* initialize a0 (seed)  		*/
	r0 = ALO		/* points to alo			*/
	rb = ALO		/* beginning of circular memory buffer 	*/
	re = CLO		/* end of circular memory buffer	*/
	pt = alo		/* point to table of constants		*/
	do 3 {
	    y=a0 x=*pt++	/* save constants alo, chi, and clo to	*/
		*r0++ = x	/* memory locations ALO, CHI, and CLO  */
	   }
	pt = alo		/* point to alo				*/
	y=a0 x=*pt++i		/* x will always contain constant alo	*/

/* NOW WE ARE SET UP, AND READY TO ROLL OUT THE RANDOM NUMBERS QUICKLY */

do n {	/* Execute the loop in the instruction cache (needed for speed) */
	y = a0			/* y = Rhi				*/
	a1l = *r0++		/* a1l = alo (zz for Rlo >= 32768)	*/
	a0 = a0 << 16		/* a0h = Rlo				*/
	if pl a1=a1<<16		/* clear zz if Rlo < 32768		*/
	p=x*y y=a0 x=*pt++i	/* p = alo * Rhi; y = Rlo		*/
	a1 = a1 + p		/* 2nd partial product + zz		*/
	a1 = a1 << 16		/* save only 16 LS bits in a1h		*/
	a0 = a0 << 4		/* a0 = ahi * Rlo (3rd partial product)	*/
	p=x*y y=a0 x=*pt++i	/* p = alo*Rlo; yh = ahi*Rlo = 16*Rlo	*/
	a0=a1+y	 y =*r0++	/* add 3rd partial product, yh = cH	*/
	a0=a0+p	 yl=*r0++	/* add 1st partial product, yl = cL	*/
	a0=a0+y			/* add c				*/
	pdx0 = a1		/* output high part of R		*/
   }
__________________________________________________________________________

For the PIC:
In the example below it computes: Rnew = (Rold * 221 + 53) mod 256

;221 = 256 - 32 - 4 + 1
;256 can be eliminated
;so we need to calculate Rnew = Rold * (1 - 32 - 4) + 53 using
;truncating arithmetic
;or Rnew = Rold * (-32 - 3) + 53

        clrc
        rlf RANDOM, f
        swapf RANDOM, w
        andlw 0xE0
        rrf RANDOM, f
        addwf RANDOM, w
        addwf RANDOM, w
        addwf RANDOM, w
        sublw 0x35
        movwf RANDOM

  .